Basic concept of programming, data types, vectors and factors
Matrices and data.frames, data import/export, graphics and visualization
Tidy data, data wrangling
Today: Spatial data in R
Basic concept of programming, data types, vectors and factors
Matrices and data.frames, data import/export, graphics and visualization
Tidy data, data wrangling
Today: Spatial data in R
library(raster)
ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100) ras
## class : RasterLayer ## dimensions : 100, 100, 10000 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
m <- rnorm(n = ncell(ras), mean = 0, sd = 1) values(ras) <- m plot(ras)
ras2 <- ras + 10
ras3 <- ras * ras2
ras4 <- ras^2
ras5 <- (ras - cellStats(ras, stat = mean)) /
cellStats(ras, stat = sd)
ras.na <- ras ras.na[ras.na > 0.5] <- NA
k <- matrix(c( 0, 0.5, 1,
0.5, 1.5, 2,
1.5, 6.0, 3),
ncol = 3, byrow = TRUE)
k
## [,1] [,2] [,3] ## [1,] 0.0 0.5 1 ## [2,] 0.5 1.5 2 ## [3,] 1.5 6.0 3
ras.reclass <- reclassify(x = ras, rcl = k)
ras[10, 20]
## ## -2.136716
ras[10, 20] <- 10
ras.stack <- stack(ras, ras2) ras.stack
## class : RasterStack ## dimensions : 100, 100, 10000, 2 (nrow, ncol, ncell, nlayers) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## names : layer.1, layer.2 ## min values : -3.497505, 6.502495 ## max values : 10.0000, 14.1605
layer1 <- subset(x = ras.stack, subset = 1) layer1
## class : RasterLayer ## dimensions : 100, 100, 10000 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## data source : in memory ## names : layer.1 ## values : -3.497505, 10 (min, max)
data <- as.data.frame(ras.stack) head(data)
## layer.1 layer.2 ## 1 -1.2699512 8.730049 ## 2 0.5653630 10.565363 ## 3 -0.3784131 9.621587 ## 4 -0.6302751 9.369725 ## 5 -0.1924782 9.807522 ## 6 0.2749651 10.274965
elevation <- raster("Data/elevation.envi")
elevation
## class : RasterLayer ## dimensions : 709, 1306, 925954 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : 16.90833, 27.79167, 44.33333, 50.24167 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi ## names : elevation
projection(elevation)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
See for projection definitions: http://spatialreference.org/
plot(elevation)
climate <- stack("Data/climate.envi")
nlayers(climate)
## [1] 4
names(climate)
## [1] "X.Tmax." "X.Tmin." "X.Tmean." "X.Prec."
names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")
writeRaster(elevation,
filename= "Data/elev.envi",
overwrite = TRUE)
library(sp)
coords_srs <- data.frame(x = runif(100, 17, 27),
y = runif(100, 45, 50))
pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs
## class : SpatialPoints ## features : 100 ## extent : 17.02223, 26.78149, 45.0296, 49.96563 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
proj <- projection(elevation) CRS(proj)
## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
pts_srs <- SpatialPoints(coords = coords_srs,
proj4string = CRS(proj))
plot(elevation) plot(pts_srs, add = TRUE)
library(mapview) mapview(climate)
df.climate <- extract(climate, pts_srs, df = TRUE) class(df.climate)
## [1] "data.frame"
head(df.climate)
## ID Tmax Tmin Tmean Prec ## 1 1 140.4167 32.583332 86.25000 53.66667 ## 2 2 105.2500 23.500000 64.16666 64.50000 ## 3 3 117.9167 30.000000 73.66666 68.58334 ## 4 4 145.7500 40.000000 92.75000 52.66667 ## 5 5 91.5000 9.166667 50.16667 71.58334 ## 6 6 116.3333 22.916666 69.41666 64.08334
spdf.climate <- extract(climate, pts_srs, sp = TRUE) spdf.climate
## class : SpatialPointsDataFrame ## features : 100 ## extent : 17.02223, 26.78149, 45.0296, 49.96563 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 4 ## names : Tmax, Tmin, Tmean, Prec ## min values : 61.9166679382324, -16.4166660308838, 22.5833339691162, 43.4166679382324 ## max values : 166.75, 67.9166641235352, 117, 84.5833358764648
You can acess SpatialPointsDataFrame the same way as data.frame!
head(spdf.climate)
## Tmax Tmin Tmean Prec ## 1 140.4167 32.583332 86.25000 53.66667 ## 2 105.2500 23.500000 64.16666 64.50000 ## 3 117.9167 30.000000 73.66666 68.58334 ## 4 145.7500 40.000000 92.75000 52.66667 ## 5 91.5000 9.166667 50.16667 71.58334 ## 6 116.3333 22.916666 69.41666 64.08334
spdf.climate$Tmax[1]
## [1] 140.4167
library(rgdal)
pts.sample <- readOGR('Data', 'pts_sample')
## OGR data source with driver: ESRI Shapefile ## Source: "Data", layer: "pts_sample" ## with 142 features ## It has 1 fields ## Integer64 fields read as strings: id
pts.sample
## class : SpatialPointsDataFrame ## features : 142 ## extent : 17.80209, 27.04407, 44.64419, 49.84834 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 1 ## names : id ## min values : 0 ## max values : 99
writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)